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Abstract 

We study laminar thin film fiows with large distortions in the free surface using the 
method of averaging across the flow. Two concrete problems are studied: the circular 
hydraulic jump and the flow down an inclined plane. For the circular hydraulic jump 
our method is able to handle an internal eddy and separated flow. Assuming a variable 
radial velocity profile like in Karman-Pohlhausen's method, we obtain a system of 
two ordinary differential equations for stationary states that can smoothly go through 
the jump where previous studies encountered a singularity. Solutions of the system 
are in good agreement with experiments. For the fiow down an inclined plane we 
take a similar approach and derive a simple model in which the velocity profile is not 
restricted to a parabolic or self-similar form. Two types of solutions with large surface 
distortions are found: solitary, kink-like propagating fronts, obtained when the fiow 
rate is suddenly changed, and stationary jumps, obtained, e.g., behind a sluice gate. 
We then include time-dependence in the model to study stability of these waves. This 
allows us to distinguish between sub- and supercritical fiows by calculating dispersion 
relations for wavelengths of the order of the width of the layer. 
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1 Introduction 



In this paper we develop a simple quantitative method to describe flows with a free surface 
which can undergo large distortions. Our method is capable of handling flows whose velocity 
profile may become far from parabolic — even including separation and regions of reverse 
flow. We are concerned with the case when the fluid layer is thin. For low Reynolds number 
flows the lubrication approximation can be used with great success (see e.g. |15[)- For 
high Reynolds number flows without separation an inviscid approximation and the shallow 



water equations are widely used. For moderate Reynolds numbers where these limiting 
approximations are invalid it is important to take both inertial and viscous effects into 
account in a consistent way, and yet one would like to keep the model simple enough to 
be tractable. In this paper we show that integral methods, like the ones developed by von 
Karman, can handle a class of such problems successfully. To be concrete we develop the 
method in the context of two physical examples: the circular hydraulic jump and the flow 
down an inclined plane. Both geometries support jump- or kink-like solutions with abrupt 
changes in the surface shape and internal velocity profiles. Analytical solutions for such 
flows are extremely difficult to obtain, and simple approximate theories that capture the 
phenomena are invaluable. 

The two flows are studied in separate sections, and an introduction is provided in the 
beginning of Sees. 2 and 3, respectively. In Sec. 2 we develop the theory for the circular 
hydraulic jump. We first study the boundary layer approximation to the full Navier-Stokes 
equations, and reduce it to a simple set of equations by averaging over the thickness. Sta- 
tionary solutions are obtained by solving a two-point boundary value problem for a system of 
only two ordinary differential equations. The solution is compared to previous experiments, 
showing good agreement. Taking advantage of the simplicity of the reduced equations, it 
is possible to obtain analytic approximations for the stationary solution. Two "outer" so- 
lutions connected by an "inner" transition region are studied separately and we obtain a 
relationship analogous to the shock condition in the classical shock theory, but within our 
viscous model. 

The flow down an inclined plane is then studied in Sec. 3. We use the same strategy 
as in Sec. 2 to derive a simple model for the two-dimensional flow. One family of solutions 
found in this model is kink-like traveling wave solutions that occur, e.g., when the flow rate 
is suddenly changed. Their velocity profiles along the inclined plane are found to stay close 
to parabolic even when a variable profile is assumed. There is another family of solutions 
with a sudden change in the surface that would correspond to the circular hydraulic jump in 
case of the radial geometry. These solutions can be interpreted as the stationary hydraulic 
jump, created behind a sluice gate in a river, even though turbulence is not included in the 
model. The flow downstream of the jump approaches a simple stationary flow, but the flow 
upstream is an expanding flow with a linear growth in thickness. The velocity profile departs 
considerably from parabolic near the jump. 

It is not easy to analyze the stability of the solutions with jumps obtained in Sees. 2 
and 3, even in the linear geometry. Instead, in Sec. 4, we include time-dependence in the 
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models and study the dispersion relation for the stationary flow with constant thickness. A 
well-established concept in the inviscid theory is to classify flows as super- and subcritical 
when the thickness is small and large, respectively. They do not have obvious counterparts, 
however, when viscosity is included. By looking carefully at the dispersion relation in the 
long and medium wave regime, we can classify the stationary flows into these two categories 
in our viscous model. The model shows spurious divergencies in the short wavelength region 
which we do not know how to overcome at present. This makes the model unsuited for 
direct time-dependent simulations. A short paper describing some of the main results has 
appeared earlier f^. 
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Figure 1: (a) Schematic view of the circular hydraulic jump, (b) Snapshot of a nearly 
perfect, stationary and circular hydraulic jump. Ethylene-glycol is used. 



2 The circular hydraulic jump 
2.1 Introduction to the problem 

When a jet of fluid hits a fiat horizontal surface, the fluid spreads out radially in a thin, 
rapidly flowing layer. At a certain distance from the jet a sudden thickening of the flow 
takes place, which is called the circular hydraulic jump. This is commonly seen, e.g., in the 
kitchen sink, but it is also important as a coating flow and in jet-cooling of a heated surface 
29| . In these practical flows with typically high Reynolds numbers, disturbances often make 



the jump non-stationary and distorted. In controlled laboratory experiments corresponding 
to a more moderate Reynolds number, an apparently stationary, radially-symmetric flow can 
be achieved. Such experiments was carried out by C. Ellegaard et al. and the results have 
been published elsewhere ||^, ^ |16|, We thank them for providing us with data and 

pictures. A schematic view and a video image of the circular jump are shown in Fig. |l]. 

In these experiments the hydraulic jump is formed on a flat disc with a circular rim. 
The rim height d can be varied, and is an important control parameter. Since the rim is 
located far from the impinging jet with the diameter of the disc around 36cm, it does not 
affect the jump except that it changes the height of the fluid layer hext exterior to the jump. 
The jump still forms even when (i = 0, but a larger d makes h^xt larger and, therefore, the 
jump stronger. Typically, h^xt exceeds d by l-2mm. The surface profiles for varying d are 
shown in Fig. ^. An interesting transition in the flow structure has been observed 0] as 
d is varied. For = 0, it was noticed before PP] , |T^, ^ |53[ that the jump contains an 



eddy on the bottom, called a separation bubble, whose inner edge is located very close to the 
position of the abrupt change on the surface, as illustrated in Fig. H(a). Such a hydraulic 
jump is referred to as a type / jump. While d remains small, this jump is stable, but as d 
is increased further, a wave-breaking transition occurs [§, which results in a new state of 
the flow. In this type 11 state, the flow has an additional eddy, called a roller or a surfing 
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Figure 2: Height profiles h{r) for different values of the external height hext- (The rim 
height d is controlled but not shown.) The height h{r) approaches h^xt for large values of r. 
Parameters are: the flow rate Q = 27[m£/s] and viscosity u = 7.6 x 10~^[m^/s], corresponding 
to the characteristic scales: radius r* = 2.8[cm], height = 1.4[mm], and radial velocity 
u^: = 12[cm/s]. Figure taken from 




Figure 3: A schematic picture showing two observed flow patterns: (a) type I flow, with a 
separation bubble, which occurs for small d, and (b) type II flow, with an additional roller 
eddy, for large d. Transitions between these states occur at a certain d, with surprisingly 
small hysteresis. 



wave, just under the surface as shown in Fig. 0(b). Q This state resembles a broken wave in 
the ocean, but is still apparently laminar. On reducing d, the type I pattern reappears, and 
there is almost no hysteresis associated with this transition. The transition from type I to 
II often leads also to breaking of the radial symmetry. An intriguing set of polygonal jumps 



pig , p!7|| are created rather than the circular one. In this paper we shall concentrate on the 



type I flow which already poses considerable difficulties. We hope to be able to generalize 
our approach in the future to be able to handle the transition to the type II flow. 

Considering how simple and common the circular hydraulic jump appears to be, it is 
surprising that a satisfactory systematic theory does not exist. The approach considered as 
"the standard" for the study of hydraulic jumps is to combine the inviscid shallow water 
equation with Rayleigh's shocks In the beginning of the century Lord Rayleigh treated 
P7| a discontinuity in a one-dimensional linear flow geometry. Such a structure is usually 
called a river bore if it is moving and a hydraulic jump if it is stationary and is created due 
to, e.g., variations in the river bed. His approach was based upon the analogy between the 



-'^ If d is increased even further, the jump "closes" as seen in Fig. |^. 
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shallow water theory and gas theory He assumed that, across such a shock, the mass 
and momentum flux are conserved but not the energy flux. 

In a coordinate system moving with the shock, the flow velocity vi and height hi upstream 
of the jump as well as V2 and /i2 downstream of the jump are taken to be positive constant 
values. Then, conservation of mass flux Q across the jump is given by 

Vihi = V2h2 = Q. (1) 
Conservation of momentum flux is 

hi (vl + ^ghi^ = h2 (^vl + ^gh2^ . (2) 

These shock conditions lead to the relation 

^1 2 _i+ /i + 8F| 



where Fi = yvf/ghi = {hc/hi)^^"^ is the upstream Froude number, F2 = yv2/gh2 = 
the downstream Froude number, and 

he = {Q'/gf' (4) 

is called the critical height. It is easy to see that he is always between hi and /i2, and that 
Fi > 1 > F2 if /ii < /ic < h2, and Fi < 1 < F2 ii hi > h^. > /z.2- In other words the jump 
connects a supercritical flow with F > 1 on the shallower side {h < h^.) to a subcritical flow 
with F < 1 on the deeper side {h > he). Since the Froude number measures the ratio of the 
fluid velocity v and the velocity of linear surface waves y/gh, it means that, in the moving 
frame, the flow moves more rapidly than the surface waves on the shallower side, but moves 
slower on the deeper side — in a precise analogy with the gas theory |3^. Further, it is 



found that the upstream hi must be supercritical by considering the change in the energy 
flux across the jump 



^ ^ _ 9Q{h2 - hif 

where Qe denotes the energy flux. Since the energy must be dissipated through the jump, 
i.e. Qe2 — Qei < 0, rather than generated, it is required that hi < h2. The origin of the 
dissipation is usually attributed to the turbulent motions at the discontinuity and surface 
waves carrying energy away from it. 

It is possible to apply this theory, combined with an assumption of the potential flow, 
for describing the circular hydraulic jump. However, it leads to incorrect estimates ||4^, Q 
of the radius of the jump Rj. Most notably, Rj is predicted to be sensitive to the radius of 
the impinging jet which should be greatly influenced by radius and height of the inlet nozzle 
where liquid comes out. In experiments [^, |^ such a strong tendency was not observed. 



Instead, it has been found that Rj scales with the flow rate Q with a certain power, and 
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it supports a model in which viscosity plays an important role. Watson constructed 
a model of the flow consisting of the inviscid and viscid regimes, and solved the viscid 
part assuming a similarity profile. By connecting to the specified external height hext via a 
Rayleigh shock, he obtained a prediction for the radius of the jump which compares favorably 
with the measurement |^, as we explain in Sec. 2.5. In his model the viscous layer starts 



from the stagnation point at r = on the plate and quickly reaches the surface at a small r. 
There is a fairly long stretch from this r to Rj in which the flow is fully viscous.0 Thus, one 
could neglect the inviscid region and assume a fully viscid flow everywhere in order to derive 
a simpler model. This assumption was made by Kurihara |^ and Tani |Q who started 



from the boundary layer equations developed by Prandtl 0, They took an average of 
the equations over the thickness, assuming also a similarity velocity profile. It resulted in 
a single ordinary differential equation for the stationary jump. This theory was elaborated 
in who realized that the flow outside the jump would naturally lead to a singularity at a 
large r. By identifying this singularity with the outflow over the rim of the plate, the flow 
outside the jump could be uniquely specified. By introducing a Rayleigh shock, the jump 
radius and its parameter dependence was calculated and compared to measurements. The 
model predicted the observed Rj reasonably well, as we review in Sees. 2.2-2.5. 

Obviously, treating the jump as a discontinuity provides us with no information on the 
internal structure of the jump region such as the type I to II transition of the flow patterns. 
It also seems inconsistent to assume a Rayleigh shock when viscous loss occurs in the whole 
domain. Why do we assume an extra energy loss at the "jump" where the flow is stationary 
and apparently laminar? It seems possible to attribute the energy dissipation entirely to 
laminar viscous forces, and to construct a viscous theory which produces a smooth but kink- 
like surface shape without the need for a discontinuity. Nevertheless, such a description must 
overcome a difficulty arising from the Goldstein- type singularity \W, ^ of the boundary 



layer equations in the vicinity of separation points. This singularity is thought to be an 
artificial one created by truncation of higher derivatives from the Navier-Stokes equations. 
It also arises in the "usual" boundary layer situation where a high Reynolds number flow 
passes a body, e.g., a wing. In such cases inviscid- viscid interaction is taken into account 
in order to resolve the singularity in a technique called the inverse method In our 



situation, however, there is no inviscid flow outside the layer. In Sees. 2.6-2.7 we propose a 
way to resolve the trouble in the following manner. We first include an additional degree of 
freedom in the velocity profile to make it non-self-similar, just like in the Karman-Pohlhausen 
method for the usual boundary layer theory. To describe the evolution, in r, of this free 
parameter, we couple the layer thickness to the pressure by assuming hydrostatic pressure. 
This serves as an alternative to the inverse method in the absence of a potential external 
flow. The resulting model for a stationary solution is two coupled ordinary differential 
equations, and reproduces the type I flow with a separation bubble — the one shown in 
Fig. |^(a). Comparison with the experiment is made in Sec. 2.7. It is possible to approximate 

^ This assumption is confirmed by recent laser-doppler measurements of the velocity profile before the 
jump pO| . Thus, the assumption made by JTsI , ^, ^ that the jump occurs at the point where the growing 
viscous layer touches the surface and the flow becomes fully developed, is incorrect. 
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analytically the stationary solution found in the model. In Sec. 2.8 the analysis is presented 
separately for the regions before and after the jump (i.e. two "outer" solutions) and the 
"inner" solution inside the jump region. An interesting observation on the inner solution is 
that a formal parameter /x can be introduced so that Rayleigh's shock condition is recovered 
in the limit n ^ 0. 



2.2 The full model 



We write down the complete model to describe the circular hydraulic jump under the assump- 
tion that the flow is laminar and radially symmetric without any angular velocity component. 
We take the radial and vertical coordinates f and z, and denote the velocity components by 
M and w, respectively.^ The governing equations are the continuity equation 



Uf + — + Wz 







(6) 



and the Navier-Stokes equations: 

1 
P 



Uf + UUf + WUz 



1 u 

Pf + U\ Uff + -Uf — +U-Z 



^ ^ ^ ^ ^ 1. ( ^ 1 . 
Wl + UWf + WWz = Pz — 9 + { Wff H Wf + Wzz 

p \ r 



(7) 



where subscripts denote partial differentiations such as Uf = du/di. For the boundary 
conditions we impose no-slip on the bottom: 



u{z = 0) = w{z = 0) = 0. 

The dynamic boundary conditions on the free surface z = h{t,r) are 
2pu 



P 



hjuf + W2 — 2hf (wf + Uz] 



(Jij — {wf + Uz) — 2hf 



[Uf 



Wy 



z=h 



z=h 



ak 



(8) 



(9) 







where a is the coefficient of surface tension and k is the mean local curvature of the free 
surface. We also need to satisfy the kinematic boundary condition on the free surface: 



hi + uhf 



w 



on z = h{t, r). 



(10) 



We are mostly interested in stationary solutions in this section. When the flow is stationary, 
we may integrate (0) over z from to h, and use (|1^) to obtain 



/ u{r,z)dz = q = —. 

Jo ZTT 



This quantity, the total mass flux Q or the mass flux per angle q, is a constant, given as a 
parameter in the experiment. 

We use tildes for the dimensional variables, dependent or independent. Dimensionless variables will be 
expressed by the same symbols but without tildes. In figures, however, we do not use tildes for simplicity. 
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2.3 Boundary layer approximation 

Since it is a formidable task to treat the full model as it stands, some simplifications need to 
be made. As explained in Sec. 1, the Reynolds number for the flow of the circular hydraulic 
jump is too large to justify the lubrication approximation, but is not large enough to use the 
inviscid approximation. Fortunately, the flow is "thin," i.e. runs predominantly horizontally 
along the plate. Truncation of the full model by the boundary layer approximation is quite 
natural in such a situation, and has indeed been used in previous literature p5|, |^ . In the 



boundary layer approximation pressure, viscous, and inertial terms in (|^) are all assumed to 
be of the same order, but there are only a few dominant terms in each group. For instance, 
a viscous term uuff is assumed to be negligible compared to I'Uzz- The dominant terms in 
the first equation in are determined in the usual manner: (if time-dependent), inertia 
terms uuf and wUz, the pressure term Pf/p, and the dominant viscous term uUzz- Similarly, 
from the second equation in (|^) we assume the dominant balance between pz/p and g. Here, 
unlike the usual boundary layer theory, we have taken into account the effect of gravity. 
This will couple the surface height h to the pressure, and will later turn out to be crucial for 
removing the singularities of the boundary layer approximation. 

If we denote the characteristic radius and height by and z^, respectively, then the 
second dominant balance requires the characteristic pressure to be pgz^:. Then, the first 
balance relation requires 



2 



where and are typical radial and vertical velocities, respectively, and is the charac- 
teristic time scale. The mass flux relation ( pT]) requires that 



M*r^2;^ = q. (13) 
while the continuity equation (||) requires 

^ = ^. (14) 



Solving (p!2D, (|T3D, and (|T^) uniquely determines the characteristic scales: 



{q^v-^g-^f^ ^ 2.7[cm] 



2* = {qvg ^)^^^ ~ 1.5[mm], 

= (gz/^3)i/8 _ i2[cm/s], (15) 

tf* = {q^^u^g)^^^ — 6.7[mm/s], 
U = {qjy-'g-'f^ ^ 0.22[s] 

where the estimated values correspond to a typical set of parameters used in the experiments: 
u ~ 0.1 cm^/s (for mixture of ethylene-glycol and water) and Q — 30 cm^/s, i.e., g ~ 5 
cm^/s. The values for r^. and 2;* correspond well to a typical jump radius and fluid thickness 
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in the experiments. Also, the predicted scahng can be experimentally tested by, for instance, 
measuring the dependence of the jump radius by changing the parameters such as q. In |Q 
evidence of the scaling and validity of the underlying assumption was given. 

We now use the characteristic scales (|I5|) , together with the pressure scale = pul, to 
non-dimensionalize the full equations. From (|^, we obtain 



1 u 

Ut + UUr + WUz = —Pr + Uzz + ( Mrr + -Ur ^ 



(16) 



{Wt + UWr + WWz) = —Pz — 1 + e'^Wzz + (^rr H U^r^ , 

where 

e = zjn = (q-^u^g-'Y^'. (17) 



Since e = z*/r* = 0.05 for the typical parameter values above, the assumption that the flow 
is "thin" is well satisfied, and we shall drop the terms of order and higher in the equations 
(p!6|). We also focus on stationary solutions in the rest of the section, and thus we obtain the 
simplified equations of motion: 

UUr + WUz = -Pr + Uzz /-|^gN 
= -pz-1. 

Correspondingly, within the error of O(e^), the dynamic boundary conditions (|^) are just 

(19) 



P\z=h = ^hrr 



Here we have introduced the Weber number 

W=^,=^,=f-,= ap-\q-'u'g-'fl\ (20) 

where £ = {2a / {gp)Y^'^ is the capillary length. For the parameter values above together with 
a ~ 70[dyn/cm] (maximum), we estimate that W ~ 0.01 and i ~ 3.8 [mm]. Since W is 
small, we neglect it in the study of stationary states.0'Q The second equation of ( [I8| ) and the 
first condition of (|I9|) with W set to zero yield hydrostatic pressure: 



p{r, z) = h{r) — z. (21) 



Combining ([TSD and (|2l|) , we obtain the stationary boundary layer equations: 

UUr + WUz = —h' + Uzz, (22) 

^ However, the term influences dispersion of short waves, so should be included in the stability analysis 
of stationary states, possibly together with the neglected terms of O(e^) and higher in (p^. 

^ The Reynolds number, defined as i? = u^^z^jv = {<f'v~^ g)^l^ « 18. The Reynolds number at the nozzle 
outlet is much higher, but it becomes moderate near the jump. 
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where the prime denotes the derivative with respect to r. This is supplemented by the 
dimensionless continuity equation: 

u 

Ur^ hw^ = 0, (23) 

r 

and mass flux condition: 

r / u{r,z)dz = 1. (24) 
Jo 

The boundary conditions have been reduced to: 

u{r, 0) = w{r, 0) = ^^5) 

In addition to these conditions, boundary conditions in the radial direction also need to be 
specified. We do not elaborate on them, however, since the in- and outlet conditions arise 
naturally without the need for prescription when we obtain a simplified system. 

The boundary layer equations (|2^ ) -(|25| ) form a closed system and can be solved numer- 
ically, but pose a difficulty when separated regions exist. Suppose that there is a separation 
point at r = and z = on a flat plate where the skin friction Uz vanishes. In its vicin- 
ity one finds [^] that generic solutions of (^) develop singularities of the Goldstein-type 
u ~ y/Vg — r, w ~ l/i/r^ — r. On the other hand, experiments p show separation and 
reversed fiow just behind a jump, so it is necessary to overcome this difficulty, which is well- 
known in the "usual" boundary layers around a body immersed in a high Reynolds number 
external fiow. No such singularities are observed in numerics of the full Navier-Stokes equa- 
tions in that case, and thus the trouble is thought to be due to truncation of the terms 
involving higher derivatives in r, i.e. the terms in (|16D of the order O(e^) and higher. An 
attempt to include those terms leads to intractable equations, so the inverse method ]T0|] is 
often used. In this method the feedback from the boundary layer into the external potential 
fiow is taken into account, and the coupled system is iteratively solved to remove the singu- 



larity. Without such an external flow present for the circular hydraulic jump, Higuera ^2 
has still obtained the velocity and height profiles from the boundary layer equations. His 
method, called marginal separation, is to force the boundary layer equations through the 
point of separation by choosing a special non-divergent velocity profile at the point. The 
physical reasoning for the choice of such a particular profile is rather unclear. Since our aim 
is also to obtain a simple tractable model, we have chosen a different approach. 

2.4 Averaged equations 



Rather than solving the partial differential equation ([2^ ) itself, we shall be content with 
satisfying only the mass and momentum conservation laws, derived from averaging ( P2| ) over 
the transverse 2;-direction. To do this we make an ansatz for the radial velocity profile u. One 
might expect that the singularities at separation points do not contribute to the averages 
and do not cause any harm. Such an expectation is too naive as shown in the next section. 
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since the model still shows singular behavior near the jump if the simplest velocity profile is 
assumed. Nevertheless, we show in Sec. 2.7 that the model becomes capable of going through 
the jump smoothly once enough flexibility is introduced in the assumed profile. 
We first define the average velocity at r by 



, u(r,z)dz. 
h Jo 

The total mass flux condition 
rhv = 1. 



41) can be written as 



(26) 



(27) 



Next, for each fixed r, we integrate the radial momentum equation (^) over z from to 
h{r), and use the continuity equation (^) with the surface boundary conditions (|25|) . We 
obtain the averaged momentum equation 



1 d 
rh dr 



r I u dz 





Using V and 
1 

h Jo 

we obtain 



G 



h /^\ 2 



dz. 



v{Gvy = -h' M^l^^o- 
Equations (^7p and ( PPf ) are the total mass and momentum equations. 



(28) 



(29) 



(30) 



2.5 Similarity profile for u 

The simplest assumption for the radial velocity profile is a self-similar ansatz: 
u{r, z)/v{r) = f{ri) 



(31) 



where rj = z/h{r) takes values between (bottom) and 1 (surface). Using (PB|), the 
ansatz can be rewritten in the alternate form: w{r,z) = rih'u{r,z). It is also equivalent 
to the requirement that the local inclination of the streamlines at (r, z) be proportional to 
T]h' = zh'{r)/h{r). Clearly, such an ansatz is too simple and "rigid" to describe a flow 
with separation. However, this is the assumption used in the previous literature, and we 
summarize its consequences. For more details, see 0. 
The conditions (^5]) and (p^) now imply 



m 
f(i) 



0, 



(32) 



f{r])dr] 
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They are not sufficient to uniquely determine /, and we choose one that is physically reason- 
able. Thus, a parabolic profile f{ri) = 31] — ?)/2rf' is a simple candidate. Using this choice, 
G = 6/5 is a constant from p9|), and (|30D becomes 

-vv =-h--^. (33) 

Other choices for / lead to the same equation with different numerical coefficients. Since all 
such equations, corresponding to different choices of f{i]), can be further transformed to 

vv' = -h'-^ (34) 

by suitably including numerical coefficients in the characteristic scales (0), the choice of / 
is not important in the study of qualitative behaviour and of parameter dependence. 
Using (P7D, the equation reduces to a single ordinary differential equation for v{r): 



v' [v- ^) = ^-vV. (35) 



This Kurihara-Tani equation was derived and studied in ^0[, in its dimensional form, and 
in 1^. The results can be summarized as follows. To find a solution corresponding to a 
hydraulic jump, the velocity v should be large for small r, and decrease smoothly as r 
increases. However, the model does not have such a solution. The coefficient of v' on the left 
hand side generically vanishes at some r where v' diverges. If (|35|) is solved in a parametric 
form on the (r, f)-plane, all solutions spiral around and into the fixed point {r,v) = (1, 1), 
that is a stable focus in the plane. Therefore, one must still connect solutions in the interior 
and the exterior by means of, e.g., a Rayleigh shock across which mass and momentum flux 
are conserved. When this is carried out, one finds that the shock occurs very close to r = 1 
in the dimensionless coordinates, implying that the radius of the jump in the dimensional 
coordinates scales roughly as in (|15|), i.e.: 

R, ^ {q'u-^g-^y^'. (36) 
This scaling relation (|36|) was compared to experiments |^, ^ by changing q for several 



different u. The radius of the jump indeed scaled with the mass flux g, but the exponent 
observed in the experiment was about 3/4 rather than 5/8 suggested by (|36|). To explain 
the discrepancy, Rj was calculated more accurately [Q. It was ffist proven that there is 
no solution for v{r) to the Kurihara-Tani equation that extends to r = oo. All solutions 
were found to diverge at some r = rend (constant) like h ~ {log(rcnd/''')}^^^- By identifying 
this singularity as the end of the plate where the water runs off, one may always find the 
solution of ( PSP diverging at the end of the plate of a given radius r = rgnd- By following the 
solution to smaller r, the solution before the jump and the position of the shock are uniquely 
determined assuming a connection via a Rayleigh shock. The shock location constructed in 
this way showed a good agreement P, pT| with the experiment. 
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2.6 Profile with a shape parameter 



An ansatz more flexible than (|3lD must be used for resolving the flow pattern in the vicinity 
of the jump. We shall allow the function / in (^Tj) to depend also on r. The simplest modifi- 
cation we can make is to assume / = /(r/, A(r)) so that the velocity profile is characterized by 
a single "shape parameter" A(r). The approach follows the ideas developed by von Karman 
and Pohlhausen ||3^ for the usual boundary layer flow around a body. There, separation of 
the boundary layer can occur when the pressure gradient, imposed by the external inviscid 
flow, becomes adverse. In our case, there is no external flow, but there is a pressure gradient, 
along the bottom z = 0, that is proportional to h'{r) due to the hydrostatic pressure (|2T|) . 
Thus, the possibility arises that the flow separates on 2; = near the jump where h' is large 
and pressure is increasing in r, as in the usual boundary layer flow. 

As an improvement over the parabolic profile, we approximate the velocity profile by the 
cubic: 

u{r, z) /v{r) = arj + hrf + erf , (37) 

where a, 6, c are now functions of r. Due to the boundary condition ( ]25| ) and mass flux 
condition (^7]) , the coefficients a, 6, and c can be expressed in terms of one parameter A as, 
for example: 

a = A + 3, 6 = -(5A + 3)/2, c = 4A/3. (38) 

The separation condition 

«.L=o = (39) 

is now equivalent to a = 0, or A = —3. The w-profile is parabolic when c = 0, or A = 0. 
Now that we have two unknowns h{r) and A(r), two equations are necessary. We use the 



averaged momentum equation (|30[) as the first equation. Note that G is now not a constant. 



but depends on the shape parameter A. From (^), we obtain 

G(A) = --- + —. (40) 
^ ^ 5 15 105 ^ ' 

Following the Karman-Pohlhausen choice, we choose the second equation to be the momen- 
tum equation (p2|) evaluated at ^ = 0: 



^'=«..L=o- (41) 

This connects the pressure gradient on 2; = with A. Using (^) and (|iO|), the two equations 
(0) and (^11) can be written as 

v^{G{\)v} = -h' -^{\ + ?,) 
dr ^ / ^ h^^ ' (42) 

= -1^(5^ + 3) 
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which can be simphfied to 



4A 



5X^3 (43) 

Finally, eliminating v using p7|), we obtain a nonautonomous system of two ordinary differ- 
ential equations for h{r) and A(r): 



^, ^ 5A + 3 

dG 4rX h^^{5X+^ 



(44) 



This is the model for the stationary circular hydraulic jump. It does become singular, but 
only on the lines h = and A = 7/2 which does not cause any trouble in describing a flow 
with a separated zone (A < —3). We show in the next section that the highly simplified 
model indeed contains solutions which describe the observed circular hydraulic jumps. A 
similar approach using momentum and energy conservation was used in , but they did not 
succeed in finding continuous solutions through the jump. 



2.7 Numerical solution of the integrated model 



The model (^4|) can be solved as a boundary value problem by specifying two boundary 
conditions for different values of r. Thus we impose 

(ri,/ii(ri)) and (ra, /i2(r2)), n < (45) 

where the values are taken from the measured surface height data. There is no fitting 
parameter once they are chosen, and the function h{r) and the shape parameter A(r) are 
determined. In particular, we do not need to specify the shape parameter as a part of the 
boundary conditions. This is an advantage of the simplified model since one no longer needs 
to specify the velocity profile at the inlet and/or outlet boundaries, which is not easy to do. 
In fact, we see that specifying both h and A at one r, either inside the jump or outside, and 
solving (1^) as an initial value problem is unstable. The system is extremely sensitive to the 
initial condition if one integrates ( ^f ) in the direction of increasing r from a small r or in the 
direction of decreasing r from a large r. Therefore, we choose ri and r2 near 1, typically ri 
around 0.4-0.8 and r2 around 1.2-1.6. Then, a straightforward shooting method from either 
boundary is sufficient to obtain a solution. After this is achieved, the solution is extended to 
r < ri and to r > r2 by integrating (H3|) backward from ri and forward from r2, respectively. 



Integrations in these directions are stable. 

Figure ^(a) shows two solutions of such a boundary value problem. They correspond 
to the two type I solutions in Fig. |^, reproduced here as dot-dashed curves. From each 
curve the boundary data are taken at fi = 11.8 [mm] (corresponding to dimensionless value 
ri = 0.42) and f2 = 30.0 [mm] (to r2 = 1.07). The computed solutions h{r) corresponding 
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Figure 4: (a) Two surface height profiles of type I flow, taken from the experiment in (§) 
are shown as the dot-dashed curves. Numerical solutions of the model (^) are shown as 
solid curves in both panels, and show reasonable agreement. To obtain each of the numerical 
solutions, h values were read from the experimental data at r = 11.8[mm] and r = 30.0[mm], 
then a boundary value problem was solved by the shooting method. The thick dashed 
curve represents an analytical approximation of the solutions before the jump, described in 
Sec. |2.8| .l. The formula ( ^2]) and (|53D shows good agreement with one fitting parameter, (b) 
The computed shape parameters A(r), characterizing the velocity profiles, corresponding to 
the two numerical solutions in (a). The flow is separated behind the jump where A < —3, 
and approaches the parabolic profile A = as r increases. Again, the dashed curve is an 
analytical approximation, (c) Two trajectories of (^) are shown in the {h, A)-plane. They 
correspond to solid curves in (a) and (b). 



to the data are shown in solid curves. Each curve shows a gradual decrease for small f as 
f increases, reaches a minimum at some f ^ 15 [mm], and then undergoes a sharp jump at 
f ^ 22- 23 [mm], and a slow decay after the jump. The location of the jump is about 10% off 
in each case, and the slope behind the jump is noticeably different. However, the qualitative 
behavior is well captured by the simple model. Figure ^(b) shows the shape parameter A. 
The velocity profile changes suddenly almost simultaneously with the rapid increase of the 
surface height, and a region where A < — 3, corresponding to separation, is observed in each 
case.0 The parameter A(r) recovers and appears to converge to A = (the parabolic profile) 
as r becomes large. 

The flow structure is more directly shown in Fig. |^, where the w-velocity profiles are 
computed from A at equidistant locations in r. Since magnitudes of the velocity vary a lot 
between small and large r, the profiles are scaled by the average velocity, so that the profiles 

^ If the downstream height is further reduced, however, the shape parameter A does not reach A = —3, 
and there is no separated region. Thus, our model predicts that a (weaker) jump without an eddy is possible. 
The flow near the bottom still decelerates just after the jump. 
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Figure 5: Visualization of the type I flow pattern based on the computed shape parameter 
A(r) from the model. The velocity profiles at equidistant locations in r are the horizontal 
component u, thus they are not tangential to the streamlines. Since magnitudes of the 
velocity vary greatly between small and large r, the profiles of u{r, z)/v{r) are shown. The 
streamlines separate zones which carry 10% of the flow rate. A separation bubble is present 
in the range of r where A < —3. Note the difference in the scales for the axes. The parameters 
differ from those of Fig. ^. They are: Q = 33[m£/s] and u = 1.4 x 10~^[m^/s], corresponding 
to r^, = 2.5[cm], = 1.7[mm], and = 16[cm/s]. 

of u{r, z)/v{r) are shown. The stream function ip is computed from the definition 

u = ipzh' 1 w = —iijr/r. (46) 

The dimensionless stream function varies from ■?/' = 0onz = 0to'?/' = lonz = /i. Inside 
the separated region ip < Q. The contours aX ip = —0.1, 0, 0.1 . . . , 1 are shown in the figure. 
That is, a region between two neighboring contour curves carries 10% of the mass flux. 

The surface velocity U predicted from the model is shown in Fig. The parameters are 
the ones used in Fig. |^. The model again misses the location of the jump by about 20%, 
so measurements and the curve from the model are offset, but qualitative features are well 
reproduced. The velocity outside the jump is small and decays like U oc 1/r, as can be seen 
from the log-log plot in the inset. This is consistent with an almost constant h and a nearly 
parabolic velocity profile, which we analytically demonstrate in the next section. On the 
other hand, the surface velocity decreases almost linearly before the jump. This region is 
harder to explain intuitively, but an analytical approximation is also obtained in the next 
section. At the jump a rapid, cusp-like drop in the velocity is noticed. 

Finally, we discuss the dependence of the solutions on the external height hext- Both in 
experiments and in the model the height inside the jump is little affected by the change in 
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Figure 6: Comparison of the prediction from the model with a surface velocity measurement 
by C. Ellegaard, A.E. Hansen, and A. Haaning The parameters are the same as in Fig. ^ 
Marker particles and a high-speed camera were used in order to obtain the surface velocity U 
shown as dots. The theoretical dotted curve was computed by finding a stationary solution 
h{r) and A(r) of a boundary value problem using two data points taken from the measured 
surface profile (not shown). Although the location of the jump is about 20% off, the model 
reproduces qualitative feature of the measurement very well. At small r, the velocity drops 
rapidly and almost linearly. It then shows a cusp-like drop at the jump, and decays gradually 
for large r. The final decay is proportional to 1/r as can be seen from the slope of about — 1 
in the log- log plot of the exterior region (inset). 



the external boundary condition h2{r2). The numerical solutions as well as the measured 
surface profiles in Fig. §(a,b) apparently overlap in the interior to the jump. Of course, the 
two solutions must correspond to different trajectories of the model (^) and cannot collapse 
exactly onto a single curve. However, the closeness of the solution curves in the interior to 
the jump is the cause of the difficulty of solving the initial value problem starting from a 
small r. 

If the external height is further increased, a transition from type I to II is observed in 
the experiment, as illustrated in Fig. ^ and Fig. 0. Unfortunately, no such transition is 
reproduced in the model when /12 is increased. Instead, one finds a computed solution of 
the model similar to the ones in Fig. | even for a much larger /i2- A physical mechanism 
to "break" the wave into a type II flow appears to be missing. In fact, a solution with a 
roller is prohibited by the model (^3]). The surface velocity on a roller is negative (inward). 
According to (0), the velocity at the surface is 

9 - A 

U = v{a + b + c) =v——, (47) 

where f > is the average velocity. Thus, U < iff A > 9. However, since we start with 
A ~ and the line A = 7/2 makes (^) singular, a solution with a roller is not possible. It 
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seems likely that this behavior can be traced back to the assumed pressure distribution ( pT] ) 
which does not provide any pressure gradient along the surface z = h. In a recent simulation 
of the circular hydraulic jump by Yokoi et a/.[|T| pressure buildup just behind the jump 
is observed and claimed to be crucial in breaking the jump. The non-hydrostatic pressure 
arises partly due to the surface tension in (0.1), but also due to the truncated viscous terms 
in (|1^) and (|l^). We do not know at present how best to extend our model to include the 
type II flows. 



2.8 Asymptotic analysis of the averaged system 

In this section we approximate the solutions of (|i3| ) analytically using formal perturbation 
expansions. We obtain explicit expressions for two "outer" regions: the region before the 
jump and the one after the jump. Moreover, we derive a single ordinary differential equation 
for the "inner" region near the jump. Analysis in the inner region connects a previous model 
using a Rayleigh shock with our model. 



2.8.1 Outer solution 1 (before the jump) 

First, we analyse the region before the jump where thickness of the fluid as well as the 
radius are small, compared to the exterior region. We denote the typical thickness, in the 
dimensionless coordinates, as 9, and treat it as a formal small parameter. We rescale the 
variables into H, R, and V as 

h = OH, 

r = O'^R, (48) 



and require consistent balance of the terms in ( ^41) or, equivalently, (|43|). The rescaling for v 
in the third equation of (|^), is chosen to ensure mass conservation (^Tf ) for all 6. In terms 
of the new variables, (ESl) can be written as 



(49) 



From the first equation the only consistent choice is to take a = 1/2. Then, in order to 
balance the power of 9 on both sides of the second equation, we need 

A = -3/5 + ^^Ai + . . . . (50) 

The form is also motivated by Fig. ^ in which A stays close to the value —0.6 before the 
jump. 

To find H{R) and the correction Ai, substitute ( |50D into the first equation of (|i9|). To 
the lowest order in 6 we obtain 

, / 1 rfi^ 1 \ 12 
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where G(— 0.6) = 1088/875 ~ 1.243. Solving this equation yields 
4 



R 5G'(-0.6) 



(52) 



where Ci is an arbitrary integration constant. The functional form agrees with Watson's 
self-similar solutions We also compare the lowest order term of 9 in the second equation 
of (Egl), and find that 



Ai 



RH^ dH 



5 dR 

By substituting H in (|5^) we obtain an approximate expression for A: 
'H^ 12 



A 



-1 + ^' 
5 



25G'(-0.6) 



R^H^ 



(53) 



We test the approximations ( |5^ ) and ( |55| ) in Fig. (^. The dashed curves are the theoretical 
curves of H{R) and X{R), shown in the dimensional coordinates. They match the numerical 
solutions and the measurements well before the jump. Here, the formal parameter 6 is taken 
as unity, and the one free parameter Ci was fitted to be 0.25. 



2.8.2 Outer solution 2 (after the jump) 

Let us now consider the behavior of (|43|) for large r. We again introduce a formal small 
parameter 6, but we now rescale r = O^^R. If we moreover assume that the height is of 
order 1, i.e., h = H, then the rescaling of the velocity is necessarily v = 9V due to (p^). 
Using these new variables, Eqs. (|43|) become: 



d 



4A 



dR 

dH _ _^5A + 3 



(54) 



In order to balance the terms in the first equation we choose 
A = e^\^ + .... 



(55) 



This is again consistent with the bottom panel of Fig. ^ where A apparently tends to 0, 
corresponding to the parabolic profile. Then, the terms of order unity in the second equation 
are 

dH _ 3 
Zr ~ ~ RH^ 

whose solution is 



(56) 



H 



12 log ■ 



R, 



end 



R 



1/4 



(57) 
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where -Rend is an integration constant representing the radius where the height goes to 0. 
Thus, (^31), as well as the simpler Kurihara-Tani model (^31), becomes singular when r oo. 
This seems to be a general property of models based on the boundary layer equations [Q. 
The absence of regular solutions for the system (|2^ ) - (|25|) when r — > cxd was proved in ^6 



We have attributed this lack of asymptotic solutions to the influence of the finite size of the 
plate. Indeed, a solution with vanishing height such as ( pTj ) reminds one very much of a flow 
running off the edge of a circular plate. 

The height H{R), given by equation (|57|), is a very slowly varying function of R. There is 
a long regime 1 <C i? -C -Rend where the height appears almost constant. In this intermediate 
regime the leading order of (0.1) becomes 

G(0)A f 1 ) . l^i (68) 



dRKRHj m 
where G(0) = 6/5. Therefore, 

^ - ^ - [m-^iR^ Wh) ^ IQ?\JP - ^) ■ ^^^^ 

We conclude that A(-R) oc 1/-R^ — which explains the observed approach to the parabolic 
velocity profile for large r. 

2.8.3 Inner solution near the jump: conservation of momentum 

Finally, we analyze the region around the hydraulic jump. Recall that in the Kurihara-Tani 
theory ( P5D the jump was obtained by fitting a Rayleigh shock. In this section, we show that 
our model is a natural generalization of the equation. 

To do this we return to (H), and introduce a formal parameter ^ in the left-hand side 
of the second equation. 

4{G(AM = -/,'-^(A + 3) ^^^^ 
^^h' = -(5A + 3). 

where v = l/{rh). The first equation describes the balance of inertia, hydrostatic pressure. 



and viscous forces. The value = 1 corresponds to (|42|). 

Setting /i = gives A = —0.6. Then the first equation becomes the Kurihara-Tani 
equation (|3^), except that the coefficient 6/5 = 1.2 is changed to G(— 0.6) ~ 1.243 here, since 
the profile is not parabolic. (As discussed before, the velocity profile is not so important in 
their model as long as it is self-similar.) Since our model corresponds to /x = 1, the parameter 
/i interpolates between the two models, but the correspondence of the two is not obvious 
because the limit /z = is a singular limit. We treat /x as a formal small parameter, and 
carry out a singular perturbation analysis to investigate the connection as well as to obtain 
an approximation in the jump region. 
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In Kurihara-Tani model a shock is needed to extend the solution from small to large 
values of r. Suppose the shock is situated at r = vq. Consider a small region of size /i around 
r = ro, and rescale the coordinate as r = ro + fiX. Then, in the inner coordinate X, Eq. ( |60D 
becomes 



_]__d_iG{Xl) _ _dh_ 

roh dX I r^h J dX 
5A + 3 

+ 0(/i). 



dh 
dX 



(61) 



roh^ 



We see that A = —0.6 with h an arbitrary constant are the only possible fixed points of 
(|6T|). Thus the solutions must satisfy A —0.6 for X —>■ ±oo. This correctly matches the 
external solution before the jump, but not after the jump, where A — >■ 0.[] The first equation 
can be integrated once, giving the momentum conservation. 



^(A) _^ 



(62) 



with an integration constant C3. Now we solve the second equation of (^) for A, and 

substitute it into this equation. Using (|40|) in the form G{X) = (A — 
an ordinary differential equation for h only: 



yI , we obtain 



105 



Vq/i^ dh 41' 
5 dX ^ 10 



13 r^h^ 

H h — — 

12 2 



(63) 



We look for a solution h{X) with h ^ hi as X —00 and /i ^ /i2 as X — > +00 where 
hi and /i2 are constants. Then, Eq. ( |B3D with the first boundary condition determines the 
constant C3 in terms of Tq and hi. Eliminating C3 we obtain 



105 



'roh^ dh 41' 
5 dX ^ 10 



hi 



41 
10 



h 



H (hi 

12^ 



h) - -^hih{hl - /i' 



0. 



(64) 



Plugging the second boundary condition into this equation yields a relation between hi and 
/i2, given tq. 



hihl + hlh2 -2hl = 



where 



(65) 



(66) 



Note that the singularity of the outer solution after the jump (£3)"(e1) for r — > does not allow correct 
matching for X — > +00 when ^ ^ Q. Nevertheless, our method reproduces the structure of the separation 
zone quite well. 
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is the critical height for the circular hydraulic jump.Q Solving this equation, we obtain an 
equation analogous to the shock condition (0): 

^ = i (-1 + ^l + 8{hJh)A = , ^ (67) 

It is easy to see that h^. is always between hi and h2, i.e., hi < < h2 or h2 < h^ < hi. The 
Froude number in this case could naturally be defined as F{XY = {hc/h{X)Y.f\ 

When hi is close to h^., the final height /i2 is close to h^ as well. Then, the Froude number 
is close to unity for all X, and the jump is weak, i.e. he — hi = 6 <^ 1. Then, we see from 
the balance of the terms in ( |64t ) that h = he + 6Y{6x). The leading balance reduces to 

y" = 7(1-^2) 

with 

196875 ^ 7 V/»_VS^„, 5/3 



^ = ^ll7j '■o-'Si.lrr. (68) 

Thus, in the weak jump limit, the height is given by 

h{x) = K + 5tanh(57x). (69) 

It is interesting to note that we can connect from /ii at X = — oo to /12 at X = +00 if 
hi < h2, but not if hi > /i2, just like in the Rayleigh shock. This requirement comes from 



the equation (|6^ ) self-consistently rather than making a hypothesis on the energy loss like 
we did in (|^). To see this, consider the stability of the fixed points hi and /i2 with respect to 
the governing equation ( |64D for h^ Linearizing (^) around the uniform solutions hi (where 
z = 1, 2), we obtain an equation for the perturbation 6hi in the height: 

-^(Shi) = Kidhi 
dX 

where 

2625 ro{2hl + hiihl-3hj)} 

= —1 2h!hi • 

If hi < he < h2, then i^i > > K2, showing that the fixed point h = hi is unstable and 
h = h2 stable. A trajectory departing from hi at X = —00 and arriving at /12 at X = +00 

^ In dimensional variables, the critical height is he — (G(— O.6)g^/(7fo)^^'^- This is identical to the critical 
height (^) that appeared in the Rayleigh shock, apart from the numerical factor and the influence of fo 
reflecting the radial geometry. The viscosity v only enters in the coefficient of dh/dX in the dimensional 
version of (|6^), thus does not affect he- 

^ However, it is not clear whether F defined in this way can be a measure of super- and subcriticality 
since the governing equations are not the shallow water equations and therefore propagation of disturbances 
do not obey the well-known velocity \/gh. 

^°0f course, this stability analysis is to study existence of stationary solutions, and not to study the 
stability of such solutions in the time-dependent theory. 
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Figure 7: Comparison between the full numerical solution of (^3D, the same two solutions 
as in Fig. § shown as solid curves, and solutions of the asymptotic equation (|6^) , shown as 
dashed curves. Even though the asymptotic analysis assumes —>■ 0, the solutions compare 
fairly well with the full numerics corresponding to /i = 1. The asymptotic analysis connects 
the model with the Rayleigh shock condition. See text. 



is not prohibited, and we can indeed find such a trajectory shown in Fig. |^. In contrast, if 
hi > he > h2, then the stability of the fixed points is reversed, and there is no trajectory 
going from hi to h2. 

When hi < he < h2 so that such a trajectory exists, the departure from hi is generally 
rapid, giving an impression of a "sharp corner" at the beginning of the jump, and the arrival 
at /i2 is much smoother just as shown in Fig. 0. This is because the magnitude of the stability 
coefficient Ki is large compared to that of K2- The feature is most pronounced when hi is 
small (so, /i2 is large). It vanishes as {h2 — hi) — > when Ki and K2 both tend to zero. 

In Fig. 1^ we compare solutions of (|6^) with the two solutions of the full numerical solution 
of ( ^3D shown in Fig. ^. The jump region is enlarged. Solutions of (|6^, shown as solid curves, 
are computed by fitting the values for hi and h2, and solving the equation using tq obtained 
from (p5D and (pB|). We chose an initial condition to be somewhere inside the jump, and 



integrated (0) forward and backward from it. Since (|6^) has a translational invariance with 
respect to X, the initial condition fixes the location of the jump without affecting the shapes 
of /i or A. The analysis assuming — > performs surprisingly well against the numerical 
solution for /i = 1. The size of the jump region is now of order /i, i.e., unity, and the internal 
structure is non-trivial. The single ordinary equation (p^ ) is capable of describing the eddy 
formation in this region. 
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3 Flow down an inclined plane 



3.1 Introduction to the problem 

The properties of waves running down an inclined plane is a subject of great theoretical and 
practical importance, and has attracted the attention of many researchers. Starting with 
the pioneering work of Kapitsa & Kapitsa [Q, some of the major contributions to this field 
are found in 0, |3l], ^ |ri|, [l^, ^ . The physical picture is the following. A fixed flux 



of fluid is constantly poured onto the inclined plane from above. The fluid forms a stream 
moving downwards under the action of gravity - an idealized model of a river. If the influx 
of fluid upstream is suddenly increased, it causes the height upstream to increase, and the 
extra mass of fluid to propagate downstream. In a river, this may be caused by the melting 
of snow at regions neighbouring the river's source, or by sudden rain. A river bore, on the 
other hand, is introduced at the mouth of the river by a tidal wave, for instance, and moves 
upstream. In both solitary wave can be formed, moving at a constant velocity c 

without changing its shape. 

We are particularly interested in kink-like solitary wave solutions going from one constant 
height hi to another /i2- One can identify such a solution with a heteroclinic orbit, connecting 
two stationary states The speed c depends on how much the fluid level is increased. 



i.e., the heights hi and h2. Alternatively, we can consider c as a parameter, and study the 
existence of the stationary solution h = const, depending on c. It is rather straightforward to 
see that two solutions with h = hi and h = h2 exist if c is sufficiently large. However, even if 
c is in that regime, it is hard to judge whether there exists a smooth solution connecting the 
two states. Based on the method of averaging in Sec. 2, we develop a simple model which 
helps us to derive criteria for their existence and to compute the wave form. The model also 
enables us to ask whether they appear as "Rayleigh shocks" in the sense that the flow is 
supercritical in front of the kink structure and subcritical behind it. As we shall elaborate, 
the distinction between super- and subcritical flows is a concept inherent in inviscid shallow 
water theory, and is not at all obvious for a viscous flow since now the waves will show 
dispersion as well as damping. Indeed, we find that the wave velocities corresponding to the 
largest wave lengths will always propagate both forward and backwards, as in a subcritical 
fiow. Nevertheless, if we focus on wavelengths of the order of the depth of the fiuid layer, a 
clear distinction can be made. 

There is another kind of fiow in the linear geometry in which a sudden thickening of 
height is observed. This solution is not only relevant for, e.g., the fiow of water exiting from 
a sluice but is also a direct analog of the circular hydraulic jump. The fiow streams rapidly in 
a region immediately after the sluice, and then abruptly slows down at a certain downstream 
position. It is stationary (i.e. c = 0) with a constant discharge, and is not obtained as a 
state connecting two "equilibrium" heights. In fact, the rapid fiow before the jump cannot 
be extended arbitrarily far upstream. We shall show that our models provide physically 
reasonable solutions in this case, too. 

In Sees. 3.2 and 3.3 we write down the complete system for the inclined plane problem. 
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non-dimensionalize it, simplify it using the boundary layer approximation, and average over 
the thickness in two ways. These steps are in parallel with those in Sec. 2, but we go through 
them briefly not only for completeness but also since the geometry and the characteristic 
scales are different. To seek stationary and traveling wave solutions, we write the equations 
in a coordinate frame moving at a constant speed in Sec. 3.4. Traveling waves are studied 
in detail in Sec. 3.5, and the stationary jumps in Sec. 3.6. 



3.2 The governing equations 

We consider a viscous, incompressible, two-dimensional flow. The coordinate system is 
X in the downstream direction parallel to the inclined plane, and y in the perpendicular 
direction above the plate. Denote the velocities in these directions by u{x, y, i) and w{x, y, i), 
respectively, the pressure by p{x,y,i), and the height by h{x,i). The governing equations 
for this problem are the continuity equation 

Ux + Wy = (71) 

and the Navier-Stokes equations 

. 1. . . , 

Ui + uux + wuy = — px + gsma + u [u^x + Uyy) 

\ (72) 

Wl + UWx + WWy = Py — g COS a + U {Wxx + Wyy) 

Here, a is the angle of the inclined plane (between and vr/2) measured downward from the 
horizontal line, and the subscripts denote the partial derivatives as before. The boundary 
conditions are identical to those of the radial geometry, i.e., (|8|)- ([ToD , by reading f as x and 
z as y. The local mass flux is: 

_ rh{x,i) 

q{x, t) = udy. 



Integrating the continuity equation ( [TTD in y over the thickness and using the boundary 
conditions, we obtain the flux conservation equation: 

k + qx = 0. (73) 

The equations above form a complete system apart from the inlet and outlet conditions. They 
possess a trivial stationary solution (Nusselt solution) with a constant h and the parabolic 
velocity profile: 

u[x,y,t) = ——{ri- — \, (74) 

where rj = y/h. Given this equilibrium flow, the local flux q is also uniform and steady, and 
is a function of h: 

q = udy = (75) 
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In a non-equilibrium flow we assume that the inclined plane is infinitely long, and the 
fiow sufficiently far downstream approaches this equilibrium ffow. We then treat the ffow 
rate q for x — >■ cxd as the characteristic mass ffux g*. The corresponding height h using 
([75|) is used as the length scale /i*, and = becomes the characteristic velocity. 

We non-dimensionalize the governing equations by these scales. The continuity equation is 
unchanged in form: 

Ux+Wy = 0, (76) 

and the Navier-Stokes equations become 

3 1, 

Ut + UUx + WUy = + 7^ + -7^[UxX + Uyy) 

R ^ (77) 

Wt + UWx + WWy = -Py - — h — (WXX + Wyy) 

it tana rC 

where the pressure is normalized to fml, and the Reynolds number is 



R = = — = — — ■ \<^> 



The dimensionless mass ffux is q{x,t) = hv in terms of the average velocity 
1 r>^ 

v{x, t) = T "^dy (79) 



h Jo 

whereby (p5|) becomes 

q = hv = h^ (80) 
in an equilibrium ffow of height h. 

3.3 Boundary layer equations and averaged models 

Since the ffow on the inclined plane is expected to be predominantly in the x-direction, the 
boundary layer approximation should be applicable |ll|, O] as long as separation does not 



occur. In a similar manner as the radial case, the dominant terms of ( [77D are: 



Ut + UUx + WUy = -px + ^ + ^Myy 

3 

= -Py 



3 1 



i^tana 

The dynamic boundary conditions on z = h reduce, as before, to: 

P\y=h = .... 

u\ =0 ^^^^ 
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with the Weber number in this case being 

W=-^ = — -1^. (83) 



From (pTi2) and (|8^) , the pressure is hydrostatic with contribution from the surface tension: 

p{x,y,t) = —^ {h{x,t)-y) + Wh,, (84) 

it tan a 

so, (|81|.l) becomes 

3 3 1 

Ut + UU^ + WUy = — - — + —Uyy + Wk^^^. (SS) 

tx ritana ri 



The mass conservation ( [73D is non- dimensionahzed to 

ht + {hv)^ = 0. (86) 

Now, we make an ansatz for the w-profile, and average over the thickness in order to 
obtain two simphfied models. First, we use the self-similar velocity profile: 

u{x,y,t)/v{x,t) = f{r]) (87) 

where r] = y/h{x,t) and the function /(r^) satisfies 

/(O) = 

f(i) = o 

f{r])dri = 1. 



Plug this ansatz into (|85D, multiply it by h, and average over y to obtain 

3h 3 3v 
{hv)t + G{hv\ = — - — hK - ^ + Whh,,, (89) 

together with the mass conservation (^). Here, 

G = \j\u/vfdy= ['f{v)dri 
n Jo Jo 

is a constant for a given profile in this model. We shall use G = 6/5 for concreteness, 
corresponding to the parabolic profile / = 3(?7 — rf /I). Equation ( [89| ) is the Cartesian 
analogue of the Kurihara-Tani equation (|33|), with time-dependent and surface tension terms. 

Next, we assume a variable one-parameter profile for u. As before, we use a third-order 
polynomial 

u{x,y,t) = v{x,t){ar] + br]'^ + CT]'^) (90) 

with a = A -|- 3, 6 = — (5A + 3)/2, and c = 4A/3 chosen to satisfy the conditions (|88D for 
/. The shape parameter X{x,t) is the single variable characterizing the velocity profile. To 
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describe the evolution of X{x,t) and h{x,t) we choose the same set of equations as in the 
circular hydraulic jump. The first equation is the mass flux equation (|86|) . In addition, we 
use the momentum equation (^) multiplied by h and averaged in y, and also ( pS]) evaluated 
at y = 0: 

{hv)t + {hv^G{X)), = ^ - -T^^hK - i^(A + 3) + WhK^^ 

R i?tana Rh (^91) 

= - - — K - ^(5A + 3) + WK,, 

R ft tan a Rh"^ 



where G{\) is given by (^OD as before. This system can be cast into the more compact form: 



V WR (92) 

cot a = 1 - — (5A + 3) + -—h^xx- 
Sh'^ 3 

In the following we call (H) with (H) the "similarity model"[] and (Hj) with (H) the 
"one-parameter model" . Both models inherit the trivial uniform solution from the complete 
Navier-Stokes model: h = v = q = 1, and A = (parabolic profile) for the one-parameter 
model. 



3.4 Stationary solutions in a moving coordinate frame 

Here, we are concerned with either stationary solutions or traveling waves whose surface 
profiles may show abrupt changes. Both types of solutions can be sought as stationary 
solutions in a moving coordinate system with a suitable constant velocity c, including the 
possibility c = 0. Thus, we use the traveling wave coordinate ^ = x — ct, and rewrite the 
models within this frame. 

Using the chain rule, the mass conservation (^6|) used in both models becomes 

-ch^ + {hv)^ = 

which can be integrated to 

— ch + hv = Q(const.) (93) 

where Q is the mass flux, viewed in the moving frame.0 The flow must approach the uniform 
equilibrium flow /i = 1 in the ^ ^ oo limit. Suppose it also approaches another equilibrium 
flow h = h2 in the ^ — > — oo limit. Then, using (PDD, the condition becomes 

-ch2 + hl = Q = -c+1. (94) 

Of course, /12 = 1 is a solution of this equation. In this case we might still be able to find a 
non-trivial solution of a pulse-like solitary wave form. Such solutions have previously been 

The similarity model is the "Shkadov model" considered in |l^ when W ^ 0. 

Note that the flux q{x, t) in the laboratory frame is, in general, not a constant. The discharge at the 
inlet, e.g., at a: = ~oo must be varied in time accordingly. 
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studied well [ITl], O], and we do not further seek this type of solutions. For a solution of 



other than /12 = 1, we need 

c = hl + h2 + l. (95) 
The solution that can be positive is 



-1 + 

ho = 

2 

which is positive if and only if c > 1. 

When c > 1 two different equilibrium solutions exist, and we hope to find a kink-like 
solution which connects the two limiting flows. However c > 1 is only the necessary condition 
for its existence. Sufficiency for the existence depends on the models and the parameters: R, 
a, and c. In Sec. 3.5 we shall clarify the parameter regime for finding such solutions. It turns 
out that the velocity profiles in this type of solutions do not deviate much from parabolic 
even in the one-parameter model. In this sense they correspond to somewhat "mild" jumps 
in terms of the flow structure. 

In Sec. 3.6 we find another family of solutions which approaches h = 1 as ^ ^ 00 when 
c < 1. These solutions do not start from an equilibrium state at ^ = —00. Instead, they are 
only valid for ^ larger than some value ^o- In the similarity model they are not interesting 
since they approach h = 1 smoothly. However, within the one-parameter model, an abrupt 
change is developed in both the surface and velocity profiles, sometimes with separation. We 
interpret this solution, when c = 0, as the analogue of the circular hydraulic jump in the 
Cartesian geometry. 

The presence of surface tension makes the order of the equations higher and makes it 
more difficult to compute the solutions even when they exist. We assume that W is small and 
negligible, and set W = in this section. Under this assumption we convert the averaged 
models into the moving coordinate frame at velocity c. Equation ( ^9]) in the similarity model 
becomes: 

, 6 9, 3 , , 3v 3h , . 

-c(/n,), + -M, + — = + (96) 

Using the condition (0), v can be eliminated. We obtain a first order differential equation 
for h: 

dh _ 15 {h-l){h^ + h + l-c) 

~ 'R c^h^ - 6{1 - cy + 15hy{Rtan a)' ^ ^ 

Similarly, (^) in the one-parameter model is converted to: 

^ ^ Rh (98) 

h^cota = 1 - ^(5A + 3) 

to be solved with (^). One variable, for instance v, can be eliminated so that the system 
becomes two-dimensional for h and A. 
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In the following sections we treat these averaged models as "dynamical systems", and 
view ^ as a time-like variable. Fixed points of these systems correspond to the uniform, 
equilibrium solutions of the original time-dependent equations. Note that stability in terms of 
the variable ^ is not equivalent to temporal stability of the original time-dependent equations. 



3.5 Traveling wave solutions 

Due to the relationship (|95[) which is a one-to-one map between c and /12 in the range c > 1, 
we may treat /i2 or c as the primary parameter interchangeably. Using h2 as a parameter 
corresponds physically to varying the height and discharge upstream and then observing the 
corresponding change in the wave velocity. The condition c > 1 is equivalent to /i2 > 0, and 
/i2 > 1 if c > 3. The two regimes /12 > 1 and /i2 < 1 are qualitatively different. For /12 > 1 
the discharge at ^ — — 00 is increased, and a forward-facing front travels downstream. As we 
shall see in this section, this state exists for small enough R. In contrast, /12 < 1 corresponds 
to a backward-facing front which is found to exist for large enough R but seems to us very 
likely unstable. Thus, we concentrate on the case h2 > 1 in the following.^ 



3.5.1 The similarity model 

Since (0) is a first order autonomous ordinary differential equation, the necessary condition 
for the existence of a heteroclinic orbit starting from /i2(> 1) and arriving at /i = 1 is that 
the fixed point h = 1 is stable and /i2 is unstable. By linearization, the fixed point h = 1 is 
found to be stable if 

c^-6(l-c)2 + 15/(i?tana) >0 (99) 

or, 

< 6(1 -c)^-c^ = 5hl + lOhl + U " 
where the denominator is positive for c > 3. Similarly, /i2 is found to be unstable if 

The denominator of /2 vanishes only at /i2 = h"^^^ 2.13 for the region /i2 > 1. If /i2 > h'^^^, 
then /2 < and ( |101| ) cannot be satisfied. We discard this region of /12. For 1 < /i2 < h^'^^ 
one finds that /2(/i2) > 1 > fi{h2). Thus, the necessary condition for the existence is simply 
( |10CI| ). Once the necessary condition is fulfilled, sufficiency is guaranteed. To see this, we 



If we used the geometric mean of the up- and downstream heights /i2 as the characteristic length, 
we would obtain equations whose symmetric appearance makes it easy to study the forward- and backward- 
facing fronts simultaneously. However, we have chosen to scale by the downstream height hi in order to 
treat the traveling waves as well as the stationary jumps. 
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Figure 8: Computed examples of the traveling wave solutions connecting two equilibrium 
states. Here, the angle of the plane a = 2[deg], and the height /i — /i2 = 1.5 as ^ — > — oo, 
corresponding to the front velocity c = 4.75. Three solutions for R = 3.5, 4.5, and 5.5 are 
shown, (a) Height h from solution of the similarity model (0). The front becomes steeper 
as R increases, (b) Height h from solution of the one-parameter model (0). The curves 
are quite similar to the ones in (a) except for the oscillation in the shallower side when R 
becomes close to a critical value. (See text.) (c) Shape parameter A corresponding to the 
solutions in (b). They deviate from the parabolic profile A = and oscillate (for R = 5.5), 
but only slightly. This explains the similarity between (a) and (b). 



only need to ensure that the denominator on the right hand side of (|97|) does not vanish in 
the region 1 < h < h2. Suppose it vanished at hs, then we would have 

c^/i^ - 6(1 - cY + 15/i^/(i?tana) = 0. (102) 

Comparison with (P^ ) gives us 

c^(l - hi) + 15(1 - /i^)/(i?tana) > 0. 

It is clear that hg > 1 is impossible. Thus, hg < 1, and there is no vanishing denominator 
in 1 < /i < /i2. In Fig. ^(a) we show computed solutions of (p7| ) for three different Reynolds 
numbers. The parameters a and /i2 are fixed, such that (|10CI|) becomes R < 6.95. Within 
this range, a larger R makes the propagating front sharper. 
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3.5.2 The one-parameter model 



We can eliminate v from ( p3D and (|98D , and think of trajectories on the phase portrait for 
{h, A). We look for a heteroclinic orbit starting from a fixed point {h2, 0) and arriving at (1, 0) 
as ^ ^ oo. It is necessary for its existence that the point {h2, 0) has at least one unstable 
direction and (1, 0) has at least one stable direction. Linearizing around the equilibrium 
point a.s h = he + 6h and A = + 5X, where he = 1 or h2, we obtain: 




It is straightforward to calculate the 2x2 Jacobian matrix J, and show that 

detJ= '°^'"'!;P""" (103) 

For the point {h2, 0) we have c — 3/ig = 1 + h2 — 2h\ < when h2 > I. This means that 
det J < for /i2 > 1, and the fixed point is always a saddle, having exactly one unstable 
direction. 

For the point (1, 0) we have det J > since c — Sh^ = h"^ + h2 — 2 > when /12 > 1- 
Thus, we must also compute the trace of J for h^ = 1 which can be shown to be 

tr J = + (33 - 61c + 25c^) tan a. 
R 

For the stability of (1,0) we need tr J < 0. Since 33 — 61c + 25c^ > for c > 3, this condition 
becomes 

60 60 

i^tano; < = ; -^5 —r = fJ/i2)- (104) 

33-61c + 25c2 -3- 11/12 + 14/1! + 50/ii + 25/it ^ ^ 

When this is satisfied, the fixed point is locally attracting, and a trajectory may reach it 
from any direction. Indeed, we find numerically that the condition ( [L04| ) also seems to be 



sufficient. For any R and a we have tried in the range (|104|) , a heteroclinic solution was 
found. Computed solutions for three different values of R are shown in Fig. ||(b) and (c). 
The parameters a and /12 are identical to the ones used for the similarity model in Fig. ^(a). 
The condition ( |104[ ) yields R < 5.59. The height profiles in (b) are essentially identical to 
the ones in (a). This is because the shape parameter A shown in (c) does not deviate much 
from A = 0, the parabolic profile. 

In Fig. |^(b) and (c), the solution is oscillatory around h = hi and A = for i? = 5.5. 
This is a feature seen when R becomes close to the critical value given by (|104|) . It happens 
when the type of the fixed point (1, 0) changes from a stable node to a stable focus. The 
point is a focus when det J > (trJ)^/4, which is equivalent to /+(/i2) < -Rtana < /_(/i2) 
where 

/±(^2) = , , ^ (105) 

-7 -9h2 + l6hl + 50hl + 25hi±2V5D ^ ' 
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and 



D = 2 + 3h2-9hl~ 19hl + 3hl + 15hl + 5/i^. 

It can be seen that /+(/i2) < /s(^2) < /-(^2) for /12 > 1- Therefore, a heterochnic solution 
can be found and exhibits oscillations in a small region /+(/i2) < -Rtana < /s(/i2). In 
Fig. 11(b) and (c) this condition corresponds to 4.81 < R < 5.59, so only the solution for 
i? = 5.5 shows oscillations. 



3.6 Stationary jumps 

If c < 1, the two averaged models have only one fixed point h = 1. Therefore, one might 
imagine that it is too limited to show any jump-like structures. Nevertheless, we look for 
trajectories that approach to the fixed point as ^ — > 00. Even though c = is the physically 
most interesting case, we treat the general case < c < 1. Since there is no h2, we use c as 
the prime parameter in this section. 



3.6.1 The similarity model 

The sole fixed point h = 1 must be stable to be the limiting point of a trajectory as ^ 00. 
For < c < 1, the condition is similar to ( ]99| ) but with reversed inequality 

c2-6(l-c)2 + 15/(i?tana) <0. (106) 

The singular height hg of the governing equation is still given by (|102| ), and, using a similar 
argument as before, it is easy to see that < /i^ < 1 is impossible when c < 1 . Thus, there is 
a trajectory which approaches h = 1 from below if ( |106D holds. When 1 > c > (6 — VQ)/5 ~ 
0.71, — 6(1 — c)^ > and ( |106| ) cannot be satisfied. When c < (6 — V6)/5, the condition 
is equivalent to 

15 

Rtana> — ^, (107) 

6(1 — c)^ — 

which is satisfied in a range of R tan a since the denominator of the right hand side is positive. 

Computed solutions for i? = 50, 70, and 100 are shown in Fig. ^ as dashed curves using 
a = 3[deg] and c = 0. The condition ( |106|) becomes R > 47.7, and is satisfied for all three. 
Each solution simply approaches h = 1 smoothly, clearly refiecting the first order nature 
of the model (|9^) . As ^ decreases, the height vanishes at a finite ^ and an inlet must be 
placed before this happens. If h is very small, (P^) simplifies to dh/dC, = 5/{2i?(l — c)}. The 
solution is 

for some ^ = ^0 where h = 0. There is no abrupt change in the solutions that resembles 
a stationary shock structure. If we use R smaller than the critical value, then there is no 
solution converging to h = 1. Therefore, we view the similarity model as inadequate for 
describing stationary jumps. 
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Figure 9: Computed stationary solutions for a = 3[deg] and c = 0. Dashed curves are 
solutions of the similarity model (|97|) for R = 50, 70, and 100. Solid curves are solutions of 
the one-parameter model (|98|) for R = 30, 50, and 70. A larger R corresponds to a slower 
convergence to the equilibrium flow h = 1. These solutions do not show any shock-like 
structure. 



3.6.2 The one-parameter model 

The sole fixed point of this model when c < 1 is {h,X) = (1,0). The Jacobian and its 
determinant is still given by (|103|) , but now c — 3hl = c — 3 < and, thus, det J < 0. 
Therefore, the fixed point is always a saddle in this range of c, and there is one direction 
convergent to the fixed point as ^ oo. It is easy to compute the corresponding trajectory 
by integrating backward in ^ from the vicinity of the fixed point. This solution seems to 
exist for all values of R, a and c < 1. We are interested in solutions which approach h = 1 
from below, and tend to = at some ^ = as ^ decreases. (To be physical, an inlet 
condition must be specified at some ^ > ^oO We can analyze the solutions asymptotically 
near by assuming that h ~ — ^q) as ^ — > + 0. Then, using ( ^3]) and Q = 1 — c in 
(p^), we obtain f ~ (1 — c)/{A(^ — ^o)}- Substituting these into (p8|.2) yields 

Q 43 

A ~ -0.6 + — ^(1 -Acot a){^ - ^of. 

5(1 - c) 

Finally, comparing coefficients of the dominant terms in (PSI.l) determines A as 

^_ 12 _ 1.93 

~ 5i?G'(-0.6)(l-c) ~ R{1 - c) ' 

We observe two qualitatively different types depending on the parameter values. If A 
increases at the point ^ = ^o, then the solution reaches the parabolic profile A = mono- 
tonically. This occurs when R is large, and three computed solutions are shown in Fig. ^ 
as solid curves. The height profile is qualitatively identical to the ones from the similarity 
model shown in dashed curves. They do not show any jump structure. 

On the other hand, if A decreases at S,Oi then the trajectory makes an excursion to smaller 
A, sometimes into the separation zone A < —3, before recovering toward A = 0. The condition 
to obtain the second type is Acot a > 1, or, 

i?tan« < ^^i^(l - c) ~ 1.94(1 - c) (109) 
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Figure 10: (a) Computed height h of the stationary solutions for the one-parameter model 
(^) using a = 3[deg], c = 0, and R = 5 and 10. A shock-like structure is visible, with a fast 
shooting flow in front of it and a slow equilibrium flow behind, (b) The shape parameter 
A corresponding to the solutions in (a) shows separation, A < —3, in both solutions, (c) 
Corresponding trajectories on the phase portrait of h versus A. In addition to the two 
solutions for R = 5 and 10, three more solutions for R = 20, 30, and 50 are shown. An 
excursion to small A before convergence to the fixed point at (0, 1) is visible for trajectories 
with small R. 

with G{X) given by (^Ol). Two solutions satisfying this condition are shown in Fig. p!0|(a) and 
(b). Both the height profile and the shape parameter vary in a similar manner to the one we 
obtained in the circular hydraulic jump. The phase portrait in (c) demonstrates how rapid 
and large the excursion can become for small R. This type of solution could be realized, for 
instance, as a stationary flow (c = 0) exiting a sluice gate placed at some ^ > ^o-Q 



A full-scale channel flow such as a river certainly requires a turbulence modelling, but we have been able 
to construct a miniature experimental model in which the flow remains laminar. However, our preliminary 
observation is that a pair of edge waves are created from the ends of the gate, which makes the flow three- 
dimensional. 
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4 Linear stability of equilibrium states 



It is quite difficult to carry out linear stability analysis around the stationary solutions 
and traveling wave solutions found so far. They have non-uniform profiles obtained only 
numerically and some of the solutions have singular points beyond which they cannot be 
continued. Moreover, the inlet boundary condition can strongly affect the stability properties 
of the solutions. We shall therefore focus on the linear geometry, and only study stability 
of the equilibrium flow h = const. The results are, however, expected to be applicable 
to the equilibrium flow sufficiently far downstream of the jump in the stationary solutions 
and to flows sufficiently up- and downstream of the moving front in case of the traveling 
wave solutions. Since the dispersion relation scales with the chosen characteristic length, as 
described in Sec. 4.4, we only need to consider the flow h = 1. Both the similarity model ( pOf ) 
and the one-parameter model (^) are considered, including the surface tension term which 
is expected to be relevant for stability. One of our aims is, of course, to judge when 
infinitesimal disturbances grow and whey they decay, but their propagation velocities are 
also of our great interest. By comparing the velocities to a reference velocity, which is zero 
for the stationary jump and c(> 3) for the traveling wave, we are able to classify different 
parts of the solutions as either super- or subcritical. 



4.1 Dispersion relations 

The ffist step is to linearize the models around the fixed point h = v = 1 and, for the 
one-parameter model, A = 0. We assume infinitesimal disturbances Sh, 6v, and SX, and 
decompose them into Fourier modes: 

6h,6v,6X^ e'^'^^'-'^'l (110) 

Plugging them into the linearized equations for the similarity model (^) and (|89|) , we 
obtain: 

Solving the equation the dispersion relation is found to be 
3i 6 , / — 

uj± = + -k±JDo 112 

where the discriminant is 

Do = + ?^k + 3e(^ + + Wk\ (113) 



4i?2 hR V25 i?tana. 

Similarly, from the one-parameter model ( p6|) and (p2D, we obtain the dispersion relation: 
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where 



4 . 178z, . /421 . 20 \, 2 i ,3 '20W . iRW . 



Note that this model also has only two dispersion relations, uj-^{k) and uj-{k) because the 
second equation of (|9^ ) does not include time-derivatives. 



4.2 Long wave limit 

We first study the long wave limit by taking only the lowest order terms in k. For 
the similarity model, the dispersion relation (|112D becomes 



uj+ = 3k + ik^{R - cot a) + 0{k^) 

^_ = Jl_ h - ik\R - cot a) + 0(P). ^^^^^ 
it 5 

As ^ 0, the group velocities doj^/dk 3 and duj_/dk (—3/5). Therefore, waves 
corresponding to uj_ propagate upstream, and the flow is suhcritical irrespective of R. By 
studying the dominant imaginary components of uj±, we also find that the reverse propagating 
branch uo^ is always stable, i.e. the disturbances decay, for small enough k whereas the 
forward propagating branch u;+ is stable only for small enough Reynolds number satisfying 

/?tana< 1. (117) 

The limiting dispersion is identical in the one-parameter model apart from numerical 



coefficients. For small A;, (114) becomes 



5 

a;+ = 3A; + ik'^{-R - cot a) + 0(P) 

122 A 5 (118) 

= -— - —k- ik'^i-R - cot a) + Oik^). 
5R 25 4 / V / 

Thus, the flow is always subcritical since the long waves in the u!_ branch propagate upstream 
with velocity —14/25. Again, this branch is stable for any R while the branch is stable 
only for small Reynolds numbers: 

/?tana<4/5. (119) 



4.3 Intermediate range of k 

It is quite unexpected that the flow is subcritical for any R. One would intuitively expect 
that disturbances cannot propagate upstream for sufficiently rapid flows. An explanation 
can be made by a more careful study of the dispersion relations ( |112| ) and ( |114D , or, in 
particular, the discriminants Dq and Di. 
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Figure 11: Real part of the dispersion relation showing the propagation of disturbances 
on the equilibrium flow, (a) Similarity model using (|112|) for R = 25, 30, and 35. (b) 
One-parameter model using (ITll ) for R = 20, 25, and 30. In both models a = 5[deg] 
and W = 0.01 are fixed. Three dashed and solid curves correspond to the u!+ and U- 
branches, respectively, of the dispersion relation. The u!+ has a positive slope, or group 
velocity, for all k, while the cu. branch has positive slope only when R is large. However, for 
large enough R, the region of k in which both branches have positive slopes extends from 
small k corresponding to wavelengths beyond the system size to large k with wavelengths 
smaller than the thickness of the flow. In this case the flow is essentially supercritical since 
disturbances are all carried away downstream. 



We first consider the similarity model. If the 0{k'^) term dominates in Dq, then the 
corresponding group velocities become 



duj± 




dk 5 V 25 Riana 
Both c+ and c_ become positive for 

i^tana > 5/2. 



(120) 



(121) 



We attempt to estimate such a range of k. For brevity we assume i?tana -C 25/2 so that 
the coefficient of k"^ in Dq can be approximated by 3/(i?tana). If the magnitude of the 
0{k'^) dominates in Do, then we must have 



3k' 



> 



9 



27k 



Rtana 4i?2' 



Wk^ 



that is, 



max 



/ 3 tan a 9 tan a 



AR 



< A; < 



RW tan a 



(122) 



Using R = 30, a = 5[deg], and W = 0.01, for instance, the condition ( |121| ) and (|122| ) gives 
a window 0.16 <^ <^ 10.7 in which we can hope that the O(fc^) term dominates. 

Rather than attempting a more accurate estimate of the zone, we demonstrate that such 
an interval can be in fact quite long, by plotting the real part of a;±(/c) for ( |112|) in Fig. |ll](a). 



39 



Three different values of R are used while a and W are fixed. The ujj^ branch, shown as 
dashed curves, has a positive slope for any k. Both phase and group velocities of this branch 
are positive. On the other hand, the a;_ branch, shown as solid curves, qualitatively changes 
with R. For R = 25 its slope appears to be negative for all k, indicating a subcritical flow. 
However, for a larger R there is an interval of k in which the slope becomes positive. In 
the limit k —>■ 0, the branch still has a negative slope in accordance with the analysis of the 
long wave limit in the previous section. However, the subcritical region near k = can be 
very small. One sees in Fig. |ll](a) that the curve has a positive slope already when k > 0.05 
and R = 35. The slope continues to be positive until k = 2, corresponding to a wavelength 
of half the thickness of the equilibrium flow. Since the system length is finite in practice, 
the subcritical flow in the k ^ limit cannot be achieved, and the flow becomes essentially 
supercritical for all the wave numbers observed. This defines the super- and subcritical flows 
within our viscous model, and confirms the intuitive picture of having a supercritical flow 
when the flow is sufficiently rapid. 

The situation is qualitatively identical in the one-parameter model. We obtain 



Rtana > 20/11 



and 



max 



/3tana 50 tan a 



5R 



<^ k <^ min 



60 tana 
R ' 



180 



RW tana Vi?Wtan 



a 



(123) 



(124) 



as the corresponding equations to ( |121| ) and (|122| ), respectively. Again using R = 30, 
a = 5[deg], and W = 0.01, the interval becomes 0.05 -C -C 0.18. The upper limit 
comes from the 0{k^) term in Di, and is estimated to be rather small since we have only 
compared the magnitudes. In fact, when we plot the real part of the dispersion relation 
( |114|) in Fig. |ll](b), we find that the branch has a positive group velocity for a much 
longer range of k. The supercritical flow near the k = limit is very small once again if R 
becomes as large as i? = 25. 



4.4 Super- and subcriticality for moving fronts 

The intermediate-^ behavior enables us to decide whether a given equilibrium flow is "inher- 
ently" super- or subcritical. This distinction is made based on wave velocities with respect 
to the laboratory frame. A more classical distinction of the two types arises in the context 
of the shock theory, as reviewed in Sec. 2.1. In this case velocities are measured with respect 
to a moving front; we call the flow "supercritical" if the group velocity of all the waves is 
less than the front velocity c, and "subcritical" if there is a wave component whose group 
velocity is larger than c. Here, we briefly note that the averaged equations can describe this 
traditional classification, too. 

Take a moving front such as the one shown in Fig. ^ We concentrate on the long wave 
limit k ^ 0. For ^ — >■ c» the flow approaches an equilibrium flow with h = 1. Linear waves 
propagate forward and backward with the group velocities du^/dk = 3 and duo^/dk = —3/5 



40 



according to the dispersion relation for the similarity model ( |116| ). This is a subcritical 
situation in the laboratory frame, but, since the front velocity is c = 1 + h2 + h"^ > 3, both 
these waves propagate into the front. Therefore, the flow is supercritical with respect to the 
front. 

To derive the dispersion relation of the equilibrium flow with height h2 for ,^^—00, 
consider rescaling the height by /i2- That is, we use this height as the characteristic length 
so that a wave number k must be multiplied by /i2. Since the flow rate is q2 = fro^i 
(^), the velocity has to be scaled by ^2/^2 = h^. Thus, the group velocities for this flow in 
the laboratory frame are dujj^/dk = ?)h\ and du-/dk = — (3/5)/i2- It is easy to show that 
3/i| > c = l + h2 + h\ for /12 > 1- Thus, one wave component propagates into the front while 
the other moves away from it so that the flow behind the front is subcritical. 

Therefore, the moving front has a supercritical flow on the shallower side and a subcritical 
flow on the deeper side, and can be regarded as a classical shock. Using the one-parameter 
model instead of the similarity model is qualitatively identical. 



4.5 Short wave limit 

We now come back to the stationary equilibrium flow, and study the dispersion relation in 
the short wave range. Since the derivation of the averaged equations relies on the assumption 
of predominantly horizontal flow, it is not our aim to accurately resolve wave components 
when k is large. We only hope that the short waves decay so that they do not interfere 
with meaningful dynamics when we simulate the time-dependent model. Unfortunately, the 
one-parameter model performs poorly in this respect compared to the similarity model. 



The dispersion relation of the similarity model ( |112[ ) can be approximated in the large k 
limit as 



Rew± = ±^/Wk'^ + 0{k) 



Thus, short waves in (RO) are damped out if > 0. 



If we neglect the surface tension and set = 0, the dispersion relation for large k is 
where c± is the velocity of the corresponding wave given by 



c± = - ± J- + (127) 
5 V 25 i?tana ^ ' 

Since c_ < 6/5 from ( |127| ), the branch uj_ is always stable, as can be seen from ( |12(j| ). On 
the other hand, since c+ > 6/5, the condition for the stability of the branch a;+ is c+ < 3, 
which is equivalent to 

i?tana<l. (128) 
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For a large R the equilibrium state is no longer stable, but this is reasonable in the absence 
of surface tension. 



Now, we turn into the dispersion relation of the one-parameter model ( 114 ). For large k, 
it behaves as 



cj± ~ ±k^/ytW/75 if > 0, (129) 
and as 



u± ~ ±p/yV(25tana) if = 0. (130) 

In either case one of the branches has an unstable component as ^ oo, irrespective of R 
or a. We have been unable to find a natural modification to the one-parameter model which 
prevents this unphysical behavior. Its cause may well be that the evolution of short waves 
is not well represented by the boundary layer approximation we started with. In fact, in 
the boundary layer equations (^) the higher order derivatives of x that are thought to be 
crucial for stability of the high-A; modes are neglected. In this view the similarity model (pQl) 



provides surprisingly reasonable behaviour for large k, even starting from (81). 
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5 Conclusions 



In this article we have presented a simple but fairly quantitative method of reducing flows 
with strongly deformed free surfaces to a manageable system of equations. By assuming a 
"flexible" velocity profile whose shape parameter is another dependent variable, flows with 
an internal eddy can be described. In the radial geometry our results compare well with 
experiments and we have obtained analytic expressions for the circular hydraulic jump. 

We have also studied the flow down an inclined plane. The reduced equations possess not 
only the traveling wave solutions (heteroclinic orbits) studied previously but also stationary 
jump solutions. We have found that the stationary solutions show a stronger change in the 
velocity profile than the traveling waves. 

Finally, we have classified different parts of the flows into super- and subcritical by study- 
ing the dispersion relation around the equilibrium flow. This classiflcation is standard for 
inviscid shallow water flow and in shock theory, but is is not obvious in the context of vis- 
cous flow. Indeed, for sufficiently long waves the averaged equations show that supercritical 
flow is not possible. However, waves with intermediate lengths can make the flow essentially 
supercritical. 

The only but serious defect of our reduced model which we have been unable to overcome 
is its short wavelength behavior. As it stands now, some artiflcial dissipation term to stabilize 
the short waves is necessary before time-dependent simulations are attempted. To our dismay 
a more natural treatment of this problem has so far eluded us. 
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